
\documentclass[11pt]{beamer}
\usepackage[UTF8,scheme=plain]{ctex}
\usepackage{listings}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
\usetheme{Boadilla}

\usepackage{framed} % 可以用 \begin{shaded}，即背景色块
\definecolor{shadecolor}{rgb}{0.9,0.9,0.9}

\lstset{
	basicstyle=\footnotesize,
    flexiblecolumns,
    showstringspaces = false,
	numbers=left,
    numberstyle=\tiny,
    numbersep=4pt,
	breaklines,%自动换行
	frame=shadowbox,%加边框
	backgroundcolor = \color{red!30!green!30!blue!30},
	%背景色：蓝灰色
	language=C++, % 语言为 c++
	commentstyle=\color{red!80!green!80!blue!80},% 注释用蓝灰色字体标出
	keywordstyle=\color{blue!70}, % 高亮部分用蓝色字体
	tabsize = 4,
    escapeinside=@@
}

\begin{document}
	\author{ 路毅 \hspace{0.3cm} 曲阜师范大学 }
	\date{\number\year 年 \number\month 月 \number\day 日}
	\title{计算物理第3部分：数值线性代数}
	%\subtitle{}
	%\logo{}
	%\institute{}
	%\date{}
	%\subject{}
	%\setbeamercovered{transparent}
	%\setbeamertemplate{navigation symbols}{}
	\begin{frame}[plain]
	\maketitle
\end{frame}

\kaishu

\begin{frame}
\frametitle{目录}
	\begin{itemize}
		\item {\color{blue} 第八章：线性方程组 }
		\vspace{1cm}
		\item 第九章：矩阵本征值
	\end{itemize}
\end{frame}

\begin{frame}{线性方程组}
线性方程组是一个非常常见的问题。
很多看似不相干的问题也可以转化为线性方程组求解问题。
比如常微分方程：
\begin{equation}
\frac{d^2}{dx^2} u(x) + p(x) \frac{d}{dx} u(x) +  q(x) u(x) = r(x),
\end{equation}
在$[a,b]$上成立，另有边界条件$u(a) = u(b) = 0$。

如果在$[a,b]$上取$n$个节点，将各阶导数用有限差分公式代替，
\begin{eqnarray}
\frac{d^2}{dx^2} u(x) |_{x_i} &\approx& ( f(x_{i+1}) + f(x_{i-1}) - 2f(x_i) )/h^2, \\
\frac{d}{dx} u(x)|_{x_i} &\approx& ( f(x_{i+1}) - f(x_{i-1}) )/(2h)
\end{eqnarray}
带入微分方程，就会得到$n$个线性方程，变量为$u(x_1), u(x_2), \cdots, u(x_n)$，解出来，就得到了微分方程的一个数值解。
\end{frame}

\begin{frame}{高斯消元法}
在小学三年级或者四年级，我们就学过了高斯消元法。
线性方程组如下
\begin{eqnarray}
a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n &=& b_1 \nonumber\\
a_{21} x_2 + a_{22} x_2 + \cdots + a_{2n} x_n &=& b_2 \nonumber\\
\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots &=& \cdots \nonumber\\
a_{n1} x_1 + a_{n2} x_2 + \cdots + a_{nn} x_n &=& b_n 
\end{eqnarray}
\end{frame}

\begin{frame}{高斯消元法}
也可以记作
\begin{equation}
A \vec{x} = \vec{b},
\end{equation}
其中
\begin{eqnarray}
A = \begin{pmatrix}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots & \vdots & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{pmatrix},
~~~~
\vec{x} = \begin{pmatrix}
x_1\\
x_2\\
\vdots\\
x_n
\end{pmatrix},
~~~~
\vec{b} = \begin{pmatrix}
b_1\\
b_2\\
\vdots\\
b_n
\end{pmatrix}
\end{eqnarray}
这个方程有唯一解的前提是，$A$是满秩矩阵，即$|A| \neq 0$。

\begin{itemize}
\item 如果交换A的两行，并相应交换$\vec{b}$，不影响解$\vec{x}$；
\item 如果将A的任一行乘以一个系数再加到另一行上，并对$\vec{b}$做相应操作，不影响解$\vec{x}$。
\item 如果交换A的两列，得到的解则与$\vec{x}$不一样，相当于$\vec{x}$交换相应的两行。
\end{itemize}
\end{frame}

\begin{frame}{高斯消元法}
我们可以将A与$\vec{b}$合在一块，记作$\mathcal{A}$，
\begin{equation}
\mathcal{A} \equiv \begin{pmatrix}
A & \vec{b}
\end{pmatrix}
\end{equation}
线性方程组则可以表达为
\begin{eqnarray}
\mathcal{A} \begin{pmatrix}
\vec{x} \\ -1 \end{pmatrix} = 0
\end{eqnarray}
\end{frame}

\begin{frame}{高斯消元法}

对于$\mathcal{A}$中的第$j=2,\cdots,n$行，我们做行变换，
\begin{eqnarray}
\mathcal{A}_{jk} = - \frac{ \mathcal{A}_{j1}}{ \mathcal{A}_{11}} \mathcal{A}_{1k} + \mathcal{A}_{jk}, ~~k=1,2,\cdots, n,
\label{Gauss1}
\end{eqnarray}
即用$- \frac{ \mathcal{A}_{jk}}{ \mathcal{A}_{11}}$乘以第一行，再加到第$j$行上。
这么做的后果是，第2-n行第1个元素等于零。
注意，上面的公式不是一个等式，而是一个算法操作，等号右边的$\mathcal{A}$是操作前的矩阵，等号左边的$\mathcal{A}$是更新后的矩阵。

类似地，我们可以将第3-n行第2个元素消为零，$\forall j = 2,3,\cdots,n$，
\begin{eqnarray}
\mathcal{A}_{jk} = - \frac{ \mathcal{A}_{j2}}{ \mathcal{A}_{22}} \mathcal{A}_{2k} + \mathcal{A}_{jk}, ~~k=2,\cdots, n,
\label{Gauss2}
\end{eqnarray}
\end{frame}

\begin{frame}{高斯消元法}

如此这般，进行下去，第$i$组行变换为\\
\vspace{0.5cm}
$\forall j = i,i+1,\cdots,n$，
\begin{eqnarray}
\mathcal{A}_{jk} = - \frac{ \mathcal{A}_{j,i}}{ \mathcal{A}_{i,i}} \mathcal{A}_{i,k} + \mathcal{A}_{j,k}, ~~k=i,\cdots, n,
\label{Gaussi}
\end{eqnarray}
最终$\mathcal{A}$会变成一个上三角矩阵：
\begin{equation}
\mathcal{A} = \begin{pmatrix}
\mathcal{A}_{11} & \mathcal{A}_{12} &  \cdots & \mathcal{A}_{1,n} & \mathcal{A}_{1,n+1} \\
0                & \mathcal{A}_{22} &  \cdots & \mathcal{A}_{2,n} & \mathcal{A}_{2,n+1} \\
0                & 0                &  \cdots & \mathcal{A}_{3,n} & \mathcal{A}_{3,n+1} \\
\vdots           & \vdots           &  \vdots & \vdots            & \vdots \\
0                & 0                &  \cdots & \mathcal{A}_{n,n} & \mathcal{A}_{n,n+1} 
\end{pmatrix}
\end{equation}
\end{frame}

\begin{frame}{高斯消元法}
变换到这一步，求解 $\mathcal{A} \begin{pmatrix}
\vec{x} \\ -1
\end{pmatrix}
=0$ 就会变得非常简单：$\forall j = n, n-1, \cdots, 1$，
\begin{eqnarray}
x_j = \frac{1}{ \mathcal{A}_{j,j} } \{ \mathcal{A}_{j,n+1} - \sum^n_{k=j+1} \mathcal{A}_{jk} x_k \}
\end{eqnarray}
\end{frame}

\begin{frame}{高斯主元消元法}
但是有一个地方可能出现问题，比如（\ref{Gauss1}）或（\ref{Gauss2}）中，如果$\mathcal{A}_{11}$或$\mathcal{A}_{22}$等于零，或者是一个非常小的双精度数，就无法操作，或者带来很大的误差。
所以我们需要在第k组行变换之前，增加一步：
\begin{shaded*}
对第$k,\cdots,n$行进行行变换之前，如果$\mathcal{A}_{k,k}$的绝对值太小，就寻找$\mathcal{A}_{k,k}, \mathcal{A}_{k+1,k}, \cdots, \mathcal{A}_{n,k}$中的最大值，并交换相应的行与第$k$行。
\end{shaded*}
这就叫做寻找{\bf 主元}，相应地，这样做叫做高斯主元消去法。

如果我们把加减乘除看做是一次“操作”，高斯消元法解一个n元线性方程组需要$O(n^3)$次操作。
所以，如果$n=1000$，电脑主频在$GHz$数量级，解这样的线性方程组可能需要几秒钟。
\end{frame}

\begin{frame}[fragile]{高斯主元消元法}
下面的代码实现高斯主元消去法。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>

/*
 * GaussLinear: 求解n元线性方程组
 * 求出的解储存在x中
 * 如果主元的绝对值小于min，就会触发行交换
 * 如果A的行列式为零，行交换会失败，程序会报错并异常退出
 */
void GaussLinear(int n, double **A, double *x, double min){

        int i,j,k;
        double y;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{高斯主元消元法}
\lstset{firstnumber=last}
\begin{lstlisting}
        for(i=0;i<n;i++){
                if( fabs( A[i][i] ) < min ){// 主元太小，触发行交换
                        for(j=i+1;j<n;j++){
                                if( fabs( A[j][i] ) > min ){
                                        for(k=i;k<n;k++){//交换第j行与第i行
                                                y = A[i][k];
                                                A[i][k] = A[j][k];
                                                A[j][k] = y;
                                        }
                                        break;//跳出for(j=i+1...)循环
                                }
                        }
                        if(j==n){// 没有在A[i+1][i],...,A[n][i]中找到绝对值比min更大的数，说明A是奇异阵（行列式等于零）
                                cerr<<" A is singular in GaussLinear(n,A,x,min).\n";
                                exit(1);
                        }
                }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{高斯主元消元法}
\lstset{firstnumber=last}
\begin{lstlisting}
                for(j=i+1;j<n;j++){// 行变换：将第i列中第i+1,...,n个元素变为零
                        y = - A[j][i] / A[i][i];
                        for(k=0;k<n+1;k++){
                                A[j][k] = y * A[i][k] + A[j][k];
                        }
                }
                // for test only //输出第i组行变换以后的矩阵，看看是否正常，用来测试代码
                /*               //如果正常，屏幕上最终会出现上三角阵
                cout<<"i="<<i<<"\t A:"<<endl;
                for(j=0;j<n;j++){
                        for(k=0;k<n+1;k++){
                                cout<<A[j][k]<<",";
                        }
                        cout<<endl;
                }
                */
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{高斯消元法}
\lstset{firstnumber=last}
\begin{lstlisting}
        for(j=n-1;j>=0;j--){//行变换已经结束了，用上三角阵计算方程组的解
                y=0;
                for(k=j+1;k<n;k++) y += A[j][k]*x[k];
                x[j] = 1/A[j][j] * ( A[j][n] - y );
        }
}

int main(){

        int i,j,n=3;
        double B[3][4] = {{1,1,1,3},{1,2,2,5},{1,3,4,8}};
        double **A = new double * [n];//开辟A的动态内存，并赋值
        for(i=0;i<n;i++){
                A[i] = new double [n+1];
                for(j=0;j<n+1;j++){
                        A[i][j] = B[i][j];
                }
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{高斯消元法}
\begin{lstlisting}
        double *x = new double [n];
        GaussLinear(3,A,x,1E-4);//调用高斯消元法，计算方程组的解
        cout<<"x=";
        for(int i=0;i<3;i++) cout<<x[i]<<",";
        cout<<endl;
        delete [] x;//删除动态内存
        for(i=0;i<n;i++) delete [] A[i];
        delete [] A;
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{LU分解}
下三角矩阵(L矩阵)就是这样的：
\begin{equation}
L = 
\begin{pmatrix}
\begin{matrix}
\alpha_{00} & \\
\alpha_{10} & \alpha_{11}
\end{matrix}
& \text{\Large 0} \\
\begin{matrix}
\vdots & \vdots \\
\alpha_{n-1,0} & \alpha_{n-1,1}
\end{matrix}
&
\begin{matrix}
\ddots & \\
\cdots & \alpha_{n-1,n-1}
\end{matrix}
\end{pmatrix}
\end{equation}

而上三角矩阵（U矩阵）是这样的：
\begin{equation}
U = 
\begin{pmatrix}
\begin{matrix}
\beta_{00} & \beta_{01} \\
 & \beta_{11}
\end{matrix}
& \begin{matrix}
\cdots & \beta_{0,n-1} \\
\cdots & \beta_{1,n-1}
\end{matrix} \\
\text{\Large 0} 
&
\begin{matrix}
\ddots & \vdots \\
 & \beta_{n-1,n-1}
\end{matrix}
\end{pmatrix}
\end{equation}

\end{frame}

\begin{frame}[fragile]{LU分解}
任意$n\times n$的矩阵A可以分解为L矩阵与U矩阵的乘积:
\begin{eqnarray}
A &=&
\begin{pmatrix}
\begin{matrix}
a_{00} & a_{01} \\
a_{10} & a_{11}
\end{matrix}
& \begin{matrix}
\cdots & a_{0,n-1} \\
\cdots & a_{1,n-1}
\end{matrix} \\
\begin{matrix}
\vdots & \vdots \\
a_{n-1,0} & a_{n-1,1}
\end{matrix}
&
\begin{matrix}
\ddots & \vdots \\
\cdots & \alpha_{n-1,n-1}
\end{matrix}
\end{pmatrix}
\nonumber\\
&=& \begin{pmatrix}
\begin{matrix}
\alpha_{00} & \\
\alpha_{10} & \alpha_{11}
\end{matrix}
& \text{\Large 0} \\
\begin{matrix}
\vdots & \vdots \\
\alpha_{n-1,0} & \alpha_{n-1,1}
\end{matrix}
&
\begin{matrix}
\ddots & \\
\cdots & \alpha_{n-1,n-1}
\end{matrix}
\end{pmatrix}
\begin{pmatrix}
\begin{matrix}
\beta_{00} & \beta_{01} \\
 & \beta_{11}
\end{matrix}
& \begin{matrix}
\cdots & \beta_{0,n-1} \\
\cdots & \beta_{1,n-1}
\end{matrix} \\
\text{\Large 0} 
&
\begin{matrix}
\ddots & \vdots \\
 & \beta_{n-1,n-1}
\end{matrix}
\end{pmatrix} \nonumber\\
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{LU分解}
Crout算法为：
\begin{shaded*}
\begin{itemize}
\item[1] $\forall i=0,\cdots,n-1, ~~\alpha_{ii} = 1$。
\item[2] $\forall j=0,1,2,...,n-1$，即每一列，
        \begin{itemize}
        \item [i] $\forall i=0,\cdots,j, ~~\beta_{ij} = a_{ij} - \sum^{i-1}_{k=0} \alpha_{ik} \beta_{kj}$,
        \item [ii] $\forall i=j+1, \cdots, n-1,~~ \alpha_{ij} = \frac{1}{\beta_{jj}} ( a_{ij} - \sum^{j-1}_{k=0} \alpha_{ik} \beta_{kj} )$
        \end{itemize}
\end{itemize}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{LU分解举例}
\begin{equation}
A = \begin{pmatrix}
1 & 1 & 1\\
1 & 2 & 2\\
1 & 2 & 3
\end{pmatrix}
= L U = \begin{pmatrix}
\alpha_{00}\\
\alpha_{10} & \alpha_{11} \\
\alpha_{20} & \alpha_{21} & \alpha_{22}
\end{pmatrix}
\begin{pmatrix}
\beta_{00} & \beta_{01} & \beta_{02} \\
           & \beta_{11} & \beta_{12} \\
           &            & \beta_{22}\\
\end{pmatrix}
\end{equation}
Crout算法第一步：
\begin{shaded*}
\begin{itemize}
\item[1] {\color{red}$\forall i=0,\cdots,n-1, ~~\alpha_{ii} = 1$。}
\end{itemize}
\end{shaded*}
\begin{equation}
A = \begin{pmatrix}
1 & 1 & 1\\
1 & 2 & 2\\
1 & 2 & 3
\end{pmatrix}
= L U = \begin{pmatrix}
{\color{red}1}\\
\alpha_{10} & {\color{red}1} \\
\alpha_{20} & \alpha_{21} & {\color{red}1}
\end{pmatrix}
\begin{pmatrix}
\beta_{00} & \beta_{01} & \beta_{02} \\
           & \beta_{11} & \beta_{12} \\
           &            & \beta_{22}\\
\end{pmatrix}
\end{equation}
\end{frame}

\begin{frame}[fragile]{LU分解举例}
Crout算法第二步：
\begin{shaded*}
\begin{itemize}
\item[2] $\forall j={\color{red}0},1,2$，即每一列，
        \begin{itemize}
        \item [i] {\color{red}$\forall i=0,\cdots,j, ~~\beta_{ij} = a_{ij} - \sum^{i-1}_{k=0} \alpha_{ik} \beta_{kj}$},
        \item [ii] $\forall i=j+1, \cdots, n-1,~~ \alpha_{ij} = \frac{1}{\beta_{jj}} ( a_{ij} - \sum^{j-1}_{k=0} \alpha_{ik} \beta_{kj} )$
        \end{itemize}
\end{itemize}
\end{shaded*}
\begin{equation}
A = \begin{pmatrix}
1 & 1 & 1\\
1 & 2 & 2\\
1 & 2 & 3
\end{pmatrix}
= L U = \begin{pmatrix}
1\\
\alpha_{10} & 1 \\
\alpha_{20} & \alpha_{21} & 1
\end{pmatrix}
\begin{pmatrix}
{\color{red} 1} & \beta_{01} & \beta_{02} \\
           & \beta_{11} & \beta_{12} \\
           &            & \beta_{22}\\
\end{pmatrix}
\end{equation}
\end{frame}

\begin{frame}[fragile]{LU分解举例}
Crout算法第二步：
\begin{shaded*}
\begin{itemize}
\item[2] $\forall j={\color{red}0},1,2$，即每一列，
        \begin{itemize}
        \item [i] {$\forall i=0,\cdots,j, ~~\beta_{ij} = a_{ij} - \sum^{i-1}_{k=0} \alpha_{ik} \beta_{kj}$},
        \item [ii] {\color{red}$\forall i=j+1, \cdots, n-1,~~ \alpha_{ij} = \frac{1}{\beta_{jj}} ( a_{ij} - \sum^{j-1}_{k=0} \alpha_{ik} \beta_{kj} )$}
        \end{itemize}
\end{itemize}
\end{shaded*}
\begin{equation}
A = \begin{pmatrix}
1 & 1 & 1\\
1 & 2 & 2\\
1 & 2 & 3
\end{pmatrix}
= L U = \begin{pmatrix}
1\\
{\color{red} 1} & 1 \\
{\color{red} 1} & \alpha_{21} & 1
\end{pmatrix}
\begin{pmatrix}
1 & \beta_{01} & \beta_{02} \\
           & \beta_{11} & \beta_{12} \\
           &            & \beta_{22}\\
\end{pmatrix}
\end{equation}
\end{frame}

\begin{frame}[fragile]{LU分解举例}
Crout算法第二步：
\begin{shaded*}
\begin{itemize}
\item[2] $\forall j=0,{\color{red}1},2$，即每一列，
        \begin{itemize}
        \item [i] {\color{red}$\forall i=0,\cdots,j, ~~\beta_{ij} = a_{ij} - \sum^{i-1}_{k=0} \alpha_{ik} \beta_{kj}$},
        \item [ii] {\color{red}$\forall i=j+1, \cdots, n-1,~~ \alpha_{ij} = \frac{1}{\beta_{jj}} ( a_{ij} - \sum^{j-1}_{k=0} \alpha_{ik} \beta_{kj} )$}
        \end{itemize}
\end{itemize}
\end{shaded*}
\begin{equation}
A = \begin{pmatrix}
1 & 1 & 1\\
1 & 2 & 2\\
1 & 2 & 3
\end{pmatrix}
= L U = \begin{pmatrix}
1\\
{ 1} & 1 \\
{ 1} & {\color{red}1} & 1
\end{pmatrix}
\begin{pmatrix}
1          & {\color{red}1} & \beta_{02} \\
           & {\color{red}1} & \beta_{12} \\
           &            & \beta_{22}\\
\end{pmatrix}
\end{equation}
\end{frame}

\begin{frame}[fragile]{LU分解举例}
Crout算法第二步：
\begin{shaded*}
\begin{itemize}
\item[2] $\forall j=0,1,{\color{red}2}$，即每一列，
        \begin{itemize}
        \item [i] {\color{red}$\forall i=0,\cdots,j, ~~\beta_{ij} = a_{ij} - \sum^{i-1}_{k=0} \alpha_{ik} \beta_{kj}$},
        \item [ii] {\color{red}$\forall i=j+1, \cdots, n-1,~~ \alpha_{ij} = \frac{1}{\beta_{jj}} ( a_{ij} - \sum^{j-1}_{k=0} \alpha_{ik} \beta_{kj} )$}
        \end{itemize}
\end{itemize}
\end{shaded*}
\begin{equation}
A = \begin{pmatrix}
1 & 1 & 1\\
1 & 2 & 2\\
1 & 2 & 3
\end{pmatrix}
= L U = \begin{pmatrix}
1\\
{ 1} & 1 \\
{ 1} & 1 & 1
\end{pmatrix}
\begin{pmatrix}
1          & 1          & {\color{red}1} \\
           & 1          & {\color{red}1} \\
           &            & {\color{red}1}\\
\end{pmatrix}
\end{equation}
\end{frame}

\begin{frame}[fragile]{LU分解}
像上面这样做，有以下特点：
\begin{itemize}
\item [1] 每一步所需要的$\alpha, \beta$值都已经被前面的步骤计算好了。
\item [2] 每一步产生的$\alpha, \beta$值都供后面的步骤使用。
\item [3] 每个$a_{ij}$都只被调用一次，在计算$\alpha_{ij}$或$\beta_{ij}$时被调用。
\item [4] 计算第$j$列时，如果调换$A,L,U$第$j$行与第$i$行（$i>j$），不影响$L,U$矩阵前面第$0,1,\cdots,j-1$列的结果，这些结果仍然是正确的。
\end{itemize}
根据前3个特点，在代码中只需要一个$n\times n$的矩阵A，
在计算中不断用新的值覆盖旧的值，最终，A的下三角部分（不包括对角元）即$L$矩阵（其对角元均为1），上三角部分即$U$矩阵。
\end{frame}

\begin{frame}[fragile]{LU分解}
\begin{shaded*}
\begin{itemize}
\item[1] $\forall i=0,\cdots,n-1, ~~\alpha_{ii} = 1$。
\item[2] $\forall j=0,1,2,...,n-1$，即每一列，
        \begin{itemize}
        \item [i] $\forall i=0,\cdots,j, ~~\beta_{ij} = a_{ij} - \sum^{i-1}_{k=0} \alpha_{ik} \beta_{kj}$,
        \item [ii] $\forall i=j+1, \cdots, n-1,~~ \alpha_{ij} = \frac{1}{\color{red}\beta_{jj}} ( a_{ij} - \sum^{j-1}_{k=0} \alpha_{ik} \beta_{kj} )$
        \end{itemize}
\end{itemize}
\end{shaded*}
如果直接这么做，有1个漏洞，如果算法2.i中得到的$\beta_{jj}$等于0或者非常小，2.ii可能引入很大的误差。

例如：
\begin{equation}
A = \begin{pmatrix}
10^{-20} & 1 \\
1        & 1 
\end{pmatrix}
= L U 
= \begin{pmatrix}
1 & 0 \\
10^{20} & 1
\end{pmatrix}
\begin{pmatrix}
10^{-20} & 1 \\
0 & 1-10^{20}
\end{pmatrix}
\end{equation}
\end{frame}

\begin{frame}[fragile]{LU分解}

因为计算机有舍入误差，所以$1-10^{20}$会得到$-10^{20}$，
\begin{eqnarray}
A &=& \begin{pmatrix}
10^{-20} & 1 \\
1        & 1 
\end{pmatrix} = LU \\ 
&=& \begin{pmatrix}
1 & 0 \\
10^{20} & 1
\end{pmatrix}
\begin{pmatrix}
10^{-20} & 1 \\
0 & 1-10^{20}
\end{pmatrix} \\
&\Rightarrow& 
\begin{pmatrix}
1 & 0 \\
10^{20} & 1
\end{pmatrix}
\begin{pmatrix}
10^{-20} & 1 \\
0 & {\color{red}-10^{20}}
\end{pmatrix}, \\
&\Rightarrow& \begin{pmatrix}
10^{-20} & 1 \\
1    & 0
\end{pmatrix}
\neq A.
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{LU分解：稳定性的改进}
引入行变换以后，Crout算法如下：
\begin{shaded*}
\begin{itemize}
\item[1] $\forall i=0,\cdots,n-1$，计算$max(i) = |a_{ij}|_{max}, j=0,1,\cdots,n-1$。
\item[2] $\forall j=0,1,2,...,n-1$，即每一列，
        \begin{itemize}
        \item [i] $\forall i=0,\cdots,j, ~~a_{ij} = a_{ij} - \sum^{i-1}_{k=0} a_{ik} a_{kj}$,
        \item [ii] $\forall i=j+1, \cdots, n-1,~~ \alpha_{ij} =  a_{ij} - \sum^{j-1}_{k=0} a_{ik} a_{kj} $
        \item [iii] 比较\\
        $\{ |a_{jj}| / max(i), |a_{j+1,j}|/max(j+1), \cdots, |a_{n-1,j}|/max(n-1) \}$，\\
        若$i=imax$时$|\alpha_{ij}| / max(i)$最大，则交换第j行与imax行，并记下这个行变换。
        \item [iv] $\forall i=j+1, \cdots, n-1, a[i][j] /= a[j][j]$。
        \end{itemize}
\end{itemize}
\end{shaded*}
其中$max(i)$是每一行各元素绝对值的最大值，所以$\beta_{ij} / max(i)$是归一化以后的值。
做这个归一化，是因为有可能有一行所有数都非常小，我们关心的是相对的大小。
上面这些步骤做完以后，得到的并不是A的LU分解，而是A经过那些行交换（算法中已经逐一记下）以后的矩阵的LU分解。
\end{frame}

\begin{frame}[fragile]{LU分解}
做完LU分解以后，线性方程组求解就很容易了
\begin{equation}
A \vec{x} = L U \vec{x} = \vec{b},
\end{equation}
我们可以记
\begin{equation}
\vec{y} \equiv U \vec{x},
\end{equation}
则有
\begin{equation}
L \vec{y} = \vec{b},
\end{equation}
这很容易解，尤其是$L_{ii}=1,i=0,\cdots,n-1$，
\begin{shaded*}
$\forall j = 0,1,\cdots,n-1$，
\begin{equation}
y_j = b_j - \sum^{j-1}_{k=0} L_{jk} y_k,
\end{equation}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{LU分解}
解出$\vec{y}$以后，很容易根据$U\vec{x} = \vec{y}$得到$\vec{x}$，
\begin{shaded*}
$\forall j = n-1, n-2,\cdots,0$，
\begin{equation}
x_j = \frac{1}{ U_{jj} }( y_j - \sum^{n-1}_{k=j+1} U_{jk} y_k ),
\end{equation}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{LU分解}
在下面的代码中，ludcmp对系数矩阵A进行LU分解，并记下发生的行交换，ludcmp\_solve则用ludcmp的结果对$A\vec{x} = \vec{b}$进行求解。
这两个函数都参考自Numerical Recipes，做了少量改动。

\begin{lstlisting}
/*
 * ludcmp():    对非奇异方阵A进行行变换，再进行LU分解
 * 基于Crout算法，将L矩阵的所有对角元都设为1，分解过后，L矩阵的非对角元+U矩阵所有元素，都写入输入的系数矩阵double **a，所以原来的系数矩阵被破坏。如果double **a是奇异阵，程序将会报错并退出。
 * int n:           维数
 * double **a:      n x n 的系数矩阵
 * int *pivot:      记录行变换，A经过这些行变换以后变成A'，A'的第i行，实则是A的第pivot[i]行
 * int sign:        如果做了偶数次行变换，sign值为1，否则为-1
 */
void ludcmp(int n, double **a, int *pivot, int & sign){

        int i,j,k,imax;
        double big, y;
        double *max = new double [n];//每行绝对值最大的元素值
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{LU分解}
\lstset{firstnumber = last}
\begin{lstlisting}
        for(i=0;i<n;i++){
                big = 0;
                for(j=0;j<n;j++){
                        if( fabs(a[i][j]) > big ){
                                big = fabs( a[i][j] );
                        }
                }
                if( big < 1E-9 ){//如果有一行所有元素绝对值都小于1E-9，则认为double**a是奇异阵
                        cout<<"a is singular in ludcmp(n,a,pivot,sign).\n";
                        exit(1);
                }
                max[i] = big;
        }
        //cout<<"max:"; for(i=0;i<n;i++)cout<<max[i]<<","; cout<<endl;
        sign =1;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{LU分解}
\lstset{firstnumber = last,
	   backgroundcolor = \color{white}}
\begin{lstlisting}
        for(j=0;j<n;j++){//依次对L、U的第0,...,n-1列进行计算
                for(i=0;i<=j;i++){
        // 计算U矩阵第j列，计算公式为 @$\beta_{ij} = a_{ij} - \sum^{i-1}_{k=0} \alpha_{ik} \beta_{kj}$@
                        y=0;
                        for(k=0;k<= i-1;k++){
                                y+= a[i][k] * a[k][j];
                        }
                        a[i][j] -= y;
                }
                imax = j; 
                big = fabs( a[j][j] ) / max[j]; // fabs( a_{jj}) / max_j, 即 a_{ij} 在第j行中的相对大小
                for(i=j+1;i<n;i++){ // 
                        y=0;
                        for(k=0;k<=j-1;k++){
                                y+= a[i][k] * a[k][j];
                        }
                        a[i][j] -= y;
                        //计算L矩阵的第j行，计算公式为@$\alpha_{ij} = \frac{1}{\beta_{ij}}( a_{ij} - \sum^{i-1}_{k=0} \alpha_{ik} \beta_{kj} )$@，这里没有乘上@$\frac{1}{\beta_{jj}}$@
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{LU分解}
\lstset{firstnumber = last,
	   backgroundcolor = \color{white}}
\begin{lstlisting}
                        if( y = fabs( a[i][j] ) / max[i] > big ){
                                big = y;
                                imax = i;
                        }
                        //在@$ \frac{a_{j,j}}{max[j]}, \frac{a_{j+1,j}}{max[j+1]},\cdots,\frac{a_{n-1,j}}{max[n-1]}$@中, @$\frac{a_{imax,j}}{max[imax]}$@最大，所以下面交换第imax与第j行。
                }
                if( imax != j ){// 交换第j行与第imax行
                        pivot[j] = imax;
                        sign *= -1;
                        y = max[j];// max[j]与max[imax]也需要交换
                        max[j] = max[imax];
                        max[imax] = y;
                        for(k=0;k<n;k++){
                                y = a[j][k];
                                a[j][k] = a[ imax ][k];
                                a[ imax ][k] = y;
                        }
                }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{LU分解}
\lstset{firstnumber = last,
	   backgroundcolor = \color{white}}
\begin{lstlisting}
                else{
                        pivot[j] = j;
                }
                for(i=j+1;i<n;i++) a[i][j] /= a[j][j];//将L矩阵第j列的元素乘上@$\frac{1}{\beta_{jj}}$@
                //for test only 每次操作完一列，可以输出来看看，用于检查代码
                /*
                cout<<"~~~~~~~~~ j="<<j<<endl;
                for(i=0;i<n;i++){
                        for(k=0;k<n;k++){
                                cout<<a[i][k]<<",";
                        }
                        cout<<endl;
                }
                */
        }
        delete [] max;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{LU分解}
\lstset{firstnumber = last}
\begin{lstlisting}
/*
 * ludcmp_solve:  解线性方程组 Ax=b, 结果 x 写入 b
 * a 不是系数矩阵 A, 而是A经过行变换后的A'的LU分解，是ludcmp()函数的结果，行变换记录在pivot中
 * 如果U是奇异阵，程序会报错并退出。
 * int n            矩阵维数
 * double ** a      A'的LU分解，A'是A经过行变换后的矩阵
 * int *pivot:      记录行变换，A经过这些行变换以后变成A'，A'的第i行，实则是A的第pivot[i]行
 */
void ludcmp_solve(int n, double ** a, int * pivot, double * b){

        int j,k;
        double y;
        // A经过行变换变为A'，对b进行同样的行变换，可以保证方程组的根不变
        for(j=0;j<n-1;j++){
                k = pivot[j];
                y = b[j];
                b[j] = b[k];
                b[k] = y;
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{LU分解}
\lstset{firstnumber = last}
\begin{lstlisting}
        //cout<<"换行操作之后，b'= "; for(j=0;j<n;j++)cout<<b[j]<<","; cout<<endl;
        // 经过行变换以后，方程组变为A'x=b'，而A'=LU，所以我们记y=Ux，则方程组变为Ly=b'，下面的代码迭代解出y，并将y的值写入double *b，覆盖原来的值
        for(j=0;j<n;j++){
                y=b[j];
                for(k=0;k<=j-1;k++){
                        y -= a[j][k] * b[k];
                }
                b[j] = y;
        }
        //cout<<"y:"; for(j=0;j<n;j++)cout<<b[j]<<","; cout<<endl;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{LU分解}
\lstset{firstnumber = last}
\begin{lstlisting}
        // 根据 y=Ux，迭代解出x，并将x的值写入double *b，覆盖原来的值
        for(j=n-1;j>=0;j--){
                y=b[j];
                for(k=j+1;k<n;k++){
                        y -= a[j][k] * b[k];
                }
                if( fabs(a[j][j])<1E-9 ){
                        cout<<"error: U is singular in the LU decomposition of A.\n";
                        exit(1);
                }
                b[j] = y/a[j][j];
        }
        //最终double *b中的值就是方程组的解
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{行列式与矩阵}
有了LU分解，行列式马上就得到了，因为
\begin{equation}
|A| = | LU | = |L| \cdot |U|,
\end{equation}
而$|L|=1$（因为其对角元都是1），并且
\begin{equation}
|U| = \prod^{n-1}_{i=0} U_{ii},
\end{equation}
所以得到
\begin{equation}
|A| = \prod^{n-1}_{i=0} U_{ii}.
\end{equation}
如果分解的是A经过行交换以后的矩阵，那么其行列式乘以$-1/1$（奇数次/偶数次行交换），就得到A的行列式。
\end{frame}

\begin{frame}[fragile]{行列式与矩阵}
逆矩阵则可以如下求得
\begin{equation}
A A^{-1} = LU A^{-1} = E,
\end{equation}
把$A^{-1}$分为$n$个列矩阵，E也分为$n$个列矩阵，进行线性方程组求解，就得到了$A^{-1}$。
当然，在这样做之前，需要判断A的行列式是否为零。
\end{frame}

\begin{frame}[fragile]{三对角矩阵}
三对角方程组为如下形式
\begin{equation}
\begin{pmatrix}
b_1 & c_1 &    &  & \\
a_2 & b_2 & c_2 &  & \\
    & \ddots & \ddots & \ddots \\
    &     & a_{n-1} & b_{n-1} & c_{n-1} \\
    &     &       & a_n     & b_n
\end{pmatrix}
\begin{pmatrix}
u_1 \\
u_2 \\
\vdots \\
u_{n-1} \\
u_n
\end{pmatrix}
= 
\begin{pmatrix}
r_1 \\
r_2 \\
\vdots \\
r_{n-1} \\
r_n
\end{pmatrix}
\end{equation}
比如，本章一开始，遇到的常微分方程，用有限差分处理以后，得到的就是一个三对角方程组。

我们只讨论{\bf 对角占优}的三对角矩阵，这也是我们处理的问题经常拥有的属性，即
\begin{equation}
|b_j| > |a_j| + |c_j|, ~~~ j=0,\cdots,n-1
\end{equation}
前人证明了，在这个条件下，我们下面介绍的追赶法不会遇到零主元。

\end{frame}

\begin{frame}[fragile]{三对角矩阵}
如果用高斯消元法或者LU分解法，当然可以解决上面的方程组，但未免有些杀鸡用牛刀的感觉。

高斯消元法与LU分解都需要$O(n^3)$次操作。

如果用“追赶法”，则只需要$O(n)$次操作。
\begin{shaded*}
\begin{itemize}
\item [1] 令 $u_1 = r_1/b_1$，$\gamma_1 = 0$，$\forall j = 2,3,\cdots,n$，做如下运算
            \begin{itemize}
            \item [] $\gamma_j = \frac{ c_{j-1} }{ b_{j-1} - a_{j-1} \gamma_{j-1} }$
            \item [] $u_j = \frac{ r_j - a_j u_{j-1} }{ b_j - a_j \gamma_j }$
            \end{itemize}
\item [2] $\forall j = n-1, \cdots, 1$,
            \begin{itemize}
            \item [] $u_j = u_j - \gamma_{j+1} u_{j+1}$
            \end{itemize}
\end{itemize}
\end{shaded*}
第1步向前迭代，第2步向后迭代。
\end{frame}

\begin{frame}[fragile]{三对角矩阵}
下面我们证明最终结果是方程组的解。
第1步向前迭代以后，得到的$u_j$我们可以改记作为$u'_j$，则有$u'_1 = r_1/b_1$，
\begin{eqnarray}
\gamma_j &=& \frac{ c_{j-1} }{ b_{j-1} - a_{j-1} \gamma_{j-1} }
\label{tridiag1}
\\
u'_j &=& \frac{ r_j - a_j u'_{j-1} }{ b_j - a_j \gamma_j } 
\label{tridiag2}
\\
u_j &=& u'_j - \gamma_{j+1} u_{j+1}
\label{tridiag3}
\end{eqnarray}
首先，我们可以证明方程组中的第一个方程被$u_j$满足，即
\begin{equation}
b_1 u_1 + c_1 u_2 = r_1,
\end{equation}
证明如下，根据(\ref{tridiag3})，
\begin{eqnarray}
u_1 &=& u'_1 - \gamma_2 u_2 \nonumber\\
    &=& r_1/b_1 - \frac{c_1}{b_1} u_2,
\end{eqnarray}
即
\begin{equation}
b_1 u_1 + c_1 u_2 = r_1.
\end{equation}
\end{frame}

\begin{frame}[fragile]{三对角矩阵}
我们可以证明方程组中的最后一个方程，即
\begin{equation}
a_n u_{n-1} + b_n u_n = r_n.
\end{equation}
证明如下
\begin{eqnarray}
a_n u_{n-1} + b_n u_n &=& a_n ( u'_{n-1} - \gamma_n u'_n ) + b_n u'_n \nonumber\\
                    &=& a_n u'_{n-1} + (b_n - a_n \gamma_n)  \nonumber\\
                    &=& r_n,
\end{eqnarray}
上式第一个等号用了(\ref{tridiag3})，第三个等号用了(\ref{tridiag2})。
所以最后一个方程是对的。
\end{frame}

\begin{frame}[fragile]{三对角矩阵}
对于第$j=2,\cdots,n-1$个方程，将（\ref{tridiag2}）带入（\ref{tridiag3}），得到
\begin{equation}
u_j = \frac{ r_j - a_j u'_{j-1} }{ b_j - a_j \gamma_j } - \gamma_{j+1} u_{j+1},
\end{equation}
而根据（\ref{tridiag3}），
\begin{equation}
u'_{j-1} = u_{j-1} + \gamma_j u_j,
\end{equation}
将上式带入上上式，即得
\begin{equation}
(b_j - a_j \gamma_j) u_j = r_j - a_j (u_{j-1} + \gamma_j u_j) -c_j u_{j+1},
\end{equation}
整理即得
\begin{equation}
a_j u_{j-1} + b_j u_j + c_j u_{j+1} = r_j.
\end{equation}
\end{frame}

\begin{frame}[fragile]{三对角矩阵}
下面的代码实现追赶法
\begin{lstlisting}
#include<iostream>
using namespace std;

/*
 * TridiagLinear:       Solve linear equations Ax=r, where A has a tridiagonal form:
 *                      b0      c0
 *                      a1      b1      c1
 *                      0       a2      b2      c2
 *                      ...     ...     ...     ...
 *                      ...     ...     ...     ...
 *                                      a_{n-1} b_{n-1}
 * matrix a must be diagonally dominant, otherwise stability problem may appear.
 * The solution is stored in u
 */
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{三对角矩阵}
\lstset{firstnumber=last}
\begin{lstlisting}
void TridiagLinear( int n, double *a, double *b, double *c, double *u, double *r ){

        int j;
        double *gamma = new double [n];
        u[0] = r[0]/b[0];
        gamma[0] = 0;
        for(j=1;j<n;j++){
                gamma[j] = c[j-1]/(b[j-1] - a[j-1]*gamma[j-1]);
                u[j] = ( r[j] - a[j]*u[j-1] )/(b[j] - a[j]*gamma[j]);
        }
        for(j=n-2;j>=0;j--){
                u[j] -= gamma[j+1]*u[j+1];
        }
        delete [] gamma;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{三对角矩阵}
\lstset{firstnumber=last}
\begin{lstlisting}
int main(){
        int n=3;
        double *a = new double [n];
        double *b = new double [n];
        double *c = new double [n];
        double *u = new double [n];
        double *r = new double [n];
        a[0]=0; a[1]=1; a[2]=1;
        b[0]=2; b[1]=2; b[2]=2;
        c[0]=1; c[1]=1; c[2]=0;
        r[0]=3; r[1]=4; r[2]=3;

        TridiagLinear(n, a, b, c, u, r);
        cout<<"u: "; for(int i=0;i<n;i++)cout<<u[i]<<","; cout<<endl;
        return 0;
}
\end{lstlisting}
当系数矩阵$A$为对角占优，即对角元大于非对角元时，追赶法一定会得到方程组的解。
\end{frame}

\begin{frame}[fragile]{迭代法：雅克比（Jacobi）方法}
线性方程组为
\begin{equation}
A \vec{x} = \vec{b},
\end{equation}
A可以分为三部分，下三角部分（不包括对角元）、对角部分、上三角部分。
\begin{equation}
A = L + D + U,
\end{equation}
所以方程组的解应满足
\begin{equation}
(L+D+U)\vec{x} = \vec{b},
\end{equation}
即
\begin{equation}
\vec{x} = D^{-1} ( \vec{b} - (L+U)\vec{x} ),
\end{equation}
借助这个式子构造迭代式
\begin{equation}
\vec{x}^{(k+1)} = D^{-1} ( \vec{b} - (L+U)\vec{x}^{(k)} ),
\end{equation}
前人证明了，如果矩阵A是正定的，或者对角占优的，这个迭代式一定会收敛至正确解。
\end{frame}

\begin{frame}[fragile]{迭代法：高斯-席德尔（Gauss-Seidel）方法}
在雅克比方法中，迭代式为
\begin{equation}
\vec{x}^{(k+1)} = D^{-1} ( \vec{b} - (L+U)\vec{x}^{(k)} ),
\end{equation}
这里面实际上有n个迭代式，n为方程组维数。
在等号右边使用$\vec{x}$的矩阵元$x^{(k)}_i$，可能已经产生了第k+1次迭代结果$x^{(k+1)}_i$，那么，如果使用$x^{(k+1)}_i$而不是$x^{(k)}_i$，一般会加速收敛的速度。

这样做，就叫做高斯-席德尔方法。
前人证明了，如果矩阵A是正定的，或者对角占优的，这个迭代式一定会收敛至正确解。

据说，如果不满足这些严格条件，有时Jacobi和Gauss-Seidel方法也都会收敛到正确解。

迭代方法还有逐次松弛法。这里我省略了。。。
\end{frame}

\begin{frame}[fragile]{应用：最小二乘法}
现在假想有一个物理现象，所有人都很感兴趣，于是实验家、理论家都研究它。

实验家说，最有意思的可观测量是$y$，有$n$个条件都会影响它，$\vec{x}=(x_1,x_2,\cdots,x_n)~$。
所以$y$是$\vec{x}$的函数$y=f(\vec{x})$。
通过几百天夜以继日的实验，我们伟大的实验家测得了$N$组数据，即$(\vec{x}_1,y_1),~(\vec{x}_2,y_2),~\cdots,~(\vec{x}_N,y_N)~$。

理论家说，根据我最聪明的头脑，技巧高超的分析，$y$与$\vec{x}$的函数关系为
\begin{eqnarray}
y = f(\vec{x}) = \sum\limits^m_{i=1} a_i g_i(\vec{x}),
\label{hypothesis}
\end{eqnarray}
其中$g_i(\vec{x})$是已知函数，我可以用手写出来，代表着m个重要的方面。
$a_i(i=1\cdots m)~$是待定参数，可通过拟合实验数据得到。
\end{frame}

\begin{frame}[fragile]{应用：最小二乘法}
实验家对理论总是充满怀疑的，他只相信自己的数据。
但他想，就算他们拟合不好，也只会使事情更有趣。

于是他们开始做拟合。
对一组参数$~\vec{a}=(a_1,a_2,\cdots,a_m)$，$y=f(\vec{x})~$在$~N~$个数据点的误差平方和为
\begin{eqnarray}
\Delta(\vec{a})
= \sum\limits^N_{j=1} (f(\vec{x}_j)-y_j)^2
= \sum\limits^N_{j=1} (\sum\limits^m_{i=1} a_i g_i(\vec{x}_j) - y_j)^2.
\label{Delta}
\end{eqnarray}
这是一个关于$a_1, a_2, \cdots, a_m$的二次型，有且仅有一个极小值点。
\end{frame}

\begin{frame}[fragile]{应用：最小二乘法}
在极小值处，$\Delta$对任意$a_l$的一阶偏导数都等于零。
即$~\forall l \in 1,\cdots,m~$，
\begin{eqnarray}
\frac{\partial \Delta}{\partial a_l}
&=& \sum\limits^N_{j=1} 2(\sum\limits^m_{i=1} a_i g_i(\vec{x}_j) - y_j) g_l(\vec{x}_j)
\nonumber\\
&=& 2 \left\{ \sum\limits^m_{i=1} (\sum\limits^N_{j=1}g_l(\vec{x}_j) g_i(\vec{x}_j)) a_i - y_j \sum\limits^N_{j=1} g_l(\vec{x}_j) \right\} = 0 .
\end{eqnarray}
即解多元一次线性方程组
\begin{eqnarray}
\left(
\begin{aligned}
~~~~~~~~~~\cdots~~~~~~~~~~~~~~~\\
\cdots~~~
\sum\limits^N_{j=1} g_l(\vec{x}_j) g_i(\vec{x}_j)~~~\cdots\\
~~~~~~~~~~\cdots~~~~~~~~~~~~~~~
\end{aligned}
\right)
\left(
\begin{aligned}
\cdot~~\\
\cdot~~\\
a_i\\
\cdot~~\\
\cdot~~\\
\end{aligned}
\right)
=
\left(
\begin{aligned}
\cdot~~~~~~~~\\
\cdot~~~~~~~~\\
\sum\limits^N_{j=1} g_l(\vec{x}_j)y_j\\
\cdot~~~~~~~~\\
\cdot~~~~~~~~\\
\end{aligned}
\right)
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{应用：最小二乘法}
得到的解$\vec{a}$即最佳参数，将最佳参数带入$\Delta$的表达式，就得到了误差的极小值$\Delta_{min}$。

如果$\Delta_{min}$非常非常小，理论家可能会很得意。

但是如果$m$非常大，即参数个数非常多，实验家可能也不买账，毕竟，给你一百万个参数，大概什么都能拟合好。。。

所以，高贵的理论都是参数少而拟合准的。
\end{frame}

\begin{frame}[fragile]{代码：最小二乘法}
\lstset{firstnumber = last,
	   backgroundcolor = \color{white}}
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>
#include"library.h" //包括我们的程序库，就可以调用其中的 GaussLinear 函数求解线性方程组
/* LeastSquare：最小二乘法拟合
 * int N            共有N组实验数据供拟合
 * double **x_data  x_data[i] 是第i组数据的各个实验条件构成的矢量@$\vec{x}_i$@
 * double **y_data  y_data[i] 是第i组数据的@$y_i$@值
 * int n            每次实验有n个条件，即@$\vec{x}_i$@有n个分量
 * int m            理论公式有m个参数
 * double *a        理论公式参数为 a[0], ..., a[m-1]
 * double (**func)(int, double *)
 *                  理论公式为 @$f(x) = \sum^{m-1}_{i=0} a_i g_i(\vec{x})$@, func为指向函数的指针的指针，记录 @$g_0, ..., g_{m-1}$@
 * double precision: GaussLinear需要的参数，在其中，如果主元小于precision，就触发换行操作
 */
double LeastSquare(int N, double **x_data, double *y_data, int n, int m, double *a, double (**func)(int, double *), double precision ){
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{代码：最小二乘法}
\lstset{firstnumber = last,
	   backgroundcolor = \color{white}}
\begin{lstlisting}
        int l, i, j;
        double y,z;
        // double ** A_merge 是 A x = b 中系数矩阵A与b的并合矩阵，即 A_merge = (A|b)，请参考高斯消元课件，这里的 A_merge对应课件中的 @$\mathcal{A}$@
        double ** A_merge = new double * [ m ];
        for(l=0;l<m;l++){
                A_merge[l] = new double [m+1];
                // A_merge的左半部分为A，课件中解释了，@$A_{li}=\sum^{N-1}_{j=0} g_l(\vec{x}_j)g_i(\vec{x}_j)$@
                for(i=0;i<m;i++){
                        y=0;
                        for(j=0;j<N;j++){
                                y += func[l](n, x_data[j]) * func[i](n, x_data[j]);
                        }
                        A_merge[l][i] = y;
                }
                // A_merge的右半部分为b，课件中解释了，@$b_{l}=\sum^{N-1}_{j=0} g_l(\vec{x}_j)y_j$@
                y=0;
                for(j=0;j<N;j++){
                        y += func[l](n, x_data[j]) * y_data[j];
                }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{代码：最小二乘法}
\lstset{firstnumber = last,
	   backgroundcolor = \color{white}}
\begin{lstlisting}
                A_merge[l][i] = y;
        }
        // for test only 在程序调试阶段，可以手写一个简单例子，并输出A_merge，看看程序算的与手推的是否一致
        /*
        cout<<"A_merge is ready"<<endl;
        for(i=0;i<m;i++){
                for(j=0;j<m+1;j++){
                        cout<<A_merge[i][j]<<",";
                }
                cout<<endl;
        }
        */
        GaussLinear(m, A_merge, a, precision);//求解Ax=b
        //下面的代码计算@$RMSD = \sqrt{ \frac{ \sum^{N-1}_{j=0} (f(\vec{x}_j)-y_j)^2 }{N} } = \sqrt{\Delta/N}$@，其中@$\Delta$@在课件中定义了。
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{代码：最小二乘法}
\lstset{firstnumber=last}
\begin{lstlisting}
        y=0;
        for(j=0;j<N;j++){
                z=0;
                for(i=0;i<m;i++)
                        z += a[i] * func[i](n, x_data[j]);
                z -= y_data[j];
                y += z*z;
        }
        y = sqrt( y/N );
        for(i=0;i<m;i++) delete [] A_merge[i];//释放动态内存
        delete [] A_merge;
        return y;// 返还RMSD
}

double g0(int n, double *x){//返回常数
        return 1;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{代码：最小二乘法}
\lstset{firstnumber=last}
\begin{lstlisting}
double g1(int n, double *x){//返回所有x_0, ..., x_{n-1}的代数和
        double z = 0;
        for(int i=0;i<n;i++) z += x[i];
        return z;
}

double g2(int n, double *x){//返回所有x_0, ..., x_{n-1}的平方和
        double z = 0;
        for(int i=0;i<n;i++) z += x[i] * x[i];
        return z;
}

double g3(int n, double *x){//返回所有x_0, ..., x_{n-1}的三次方和
        double z = 0;
        for(int i=0;i<n;i++) z += pow(x[i],3);
        return z;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{代码：最小二乘法}
\lstset{firstnumber=last}
\begin{lstlisting}
int main(){

        int i,j,n=1, m=3,N=10;
        double ( *g[3] )(int, double * )={g0,g1,g2};//定义指向函数的指针数组，数组中元素分别指向g0, g1, g2

        double **x_data = new double * [N];// 编造“实验数据”
        for(i=0;i<N;i++) x_data[i] = new double [n];
        double *y_data = new double [N];

        for(i=0;i<N;i++){
                x_data[i][0]=i;
                y_data[i]=1.234+2.456*i+3.789*i*i;
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{代码：最小二乘法}
\lstset{firstnumber=last}
\begin{lstlisting}
        double *a = new double [m];
        
        cout<<"最优RMSD = "<<LeastSquare(N, x_data, y_data, n, m, a, g, 0.01 )<<endl;
        //调用LeastSquare函数，计算最优参数并储存在a中，并输出最优参数对应的RMSD
        cout<<"最有参数: "; for(i=0;i<m;i++)cout<<a[i]<<","; cout<<endl;

        delete [] a, y_data;
        for(i=0;i<N;i++)delete [] x_data[i]; delete [] x_data;

        return 0;
}
\end{lstlisting}                                     
\end{frame}

\begin{frame}
\frametitle{目录}
	\begin{itemize}
		\item 第八章：线性方程组
		\vspace{1cm}
		\item {\color{blue} 第九章：矩阵本征值 }
	\end{itemize}
\end{frame}

\begin{frame}[fragile]{矩阵本征值}
矩阵本征方程为
\begin{equation}
A \vec{x} = \lambda \vec{x},
\end{equation}
我们只讨论$A$为{\bf 实对称阵}的情况。\\
\vspace{0.5cm}
常用的哈密顿量是厄米的，实数的哈密顿量矩阵是实对称阵，它的本征值就是本征能量，本征函数就是体系本征态。
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
数学家们已经证明了，任意实对称阵$A$都存在正交变换$S$，使得$S^\top A S$为对角阵。

Jacobi方法构造一系列旋转，最终使得$A$近似地对角化，完成目标。
在p-q平面上逆时针旋转$\phi$角对应的变换矩阵为
\begin{equation}
P_{pq} = \begin{bmatrix}
 1 \\
  & \ddots \\
  &     & c     &     & s \\
  &     &       & \ddots \\
  &     & -s    &         & c \\
  &     &       &        &       & \ddots \\
  &     &       &        &       &       & 1
\end{bmatrix}
\end{equation}
其中$c=\cos\phi,-s,s=\sin\phi,c$分别处在第$p,q$行，第$p,q$列，矩阵的其余部分与单位矩阵相同。
注意$P_{pq}$是一个矩阵，不是带下标的矩阵元。
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
在旋转变换以后，矩阵A变为
\begin{equation}
A' = P^\top_{pq} A P_{pq}
= \begin{bmatrix}
        &        &  a'_{0p} &       & a'_{0q} &         &   \\
        &        &   \vdots &       &  \vdots &         &   \\
a'_{p0} & \cdots & a'_{pp}  & \cdots& a'_{pq} &\cdots   & a'_{p,n-1} \\
        &       & \vdots    &       & \vdots  &     \\
a'_{q0} & \cdots & a'_{qp}  & \cdots& a'_{qq} &\cdots   & a'_{q,n-1} \\  
        &        &   \vdots &       &  \vdots &         &   \\ 
        &        &  a'_{n-1,p} &       & a'_{n-1,q} &         &      
\end{bmatrix}
\end{equation}
其中，行数或者列数为p,q的矩阵元发生变化，其他矩阵元不变。
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
经过计算，矩阵元的变化如下：
\begin{eqnarray}
a'_{pp} &=& c^2 a_{pp} - 2cs a_{pq} + s^2 a_{qq},
\label{jacobi-a'pp}
\\
a'_{qq} &=& s^2 a_{pp} + 2cs a_{pq} + c^2 a_{qq},
\label{jacobi-a'qq}
\\
a'_{pq} &=& (c^2 - s^2)a_{pq} + cs( a_{pp} - a_{qq} ),
\label{jacobi-a'pq}
\\
a'_{ip} &=& c a_{ip} - s a_{iq}, ~~\forall i \neq p,q 
\label{jacobi-a'ip}
\\
a'_{iq} &=& s a_{ip} + c a_{iq}, ~~\forall i \neq p,q.
\label{jacobi-a'iq}
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
Jacobi希望让非对角元变为零，所以令$a'_{pq}=0$，即
\begin{equation}
(c^2 - s^2)a_{pq} + cs( a_{pp} - a_{qq} )=0,
\end{equation}
我们就得到了旋转角$\phi$的相关信息
\begin{eqnarray}
\theta = \cot 2 \phi = \frac{ a_{qq} - a_{pp} }{ 2 a_{pq}}.
\end{eqnarray}
如果我们记
\begin{equation}
t = \frac{s}{c} = \tan \phi,
\end{equation}
则有
\begin{equation}
t^2 + 2t\theta - 1 = 0,
\end{equation}
只要满足这个式子，旋转以后就有$a'_{pq}$。
解出来得到
\begin{equation}
t = - \theta \pm \sqrt{\theta^2 +1},
\end{equation}
\end{frame}

\begin{frame}
选取绝对值更小的t，可以使误差更小，
\begin{equation}
t = \left\{
\begin{aligned}
&&- \theta + \sqrt{\theta^2 +1} = \frac{1}{ \theta + \sqrt{\theta^2 +1}}, ~~~ \theta>0, \\
&&- \theta - \sqrt{\theta^2 +1} = - \frac{1}{ -\theta + \sqrt{\theta^2 +1}}, ~~~ \theta<0.
\end{aligned}
\right.
\end{equation}
这两个方程可以用下面的方程概括
\begin{equation}
t = \frac{ sgn(\theta) }{ |\theta| + \sqrt{\theta^2 +1}}.
\end{equation}
在这个约定下，$ |t| < 1$。
所以，我们可以从$a_{pp}, a_{pq}, a_{qq}$计算得到$\theta$，进而用上面的式子得到$t$，得到$t$以后，我们可以进一步计算
\begin{equation}
c = \frac{1}{1+t^2}, ~~ s = ct,
\end{equation}
在上式中，我们假定了$c>0$，综合$|t|<1$，我们实际上约束了$ - \frac{\pi}{4} < \phi < \frac{\pi}{4}$。
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
出于计算方便和误差控制，我们可以利用$t$，将(\ref{jacobi-a'pp}-\ref{jacobi-a'iq})改写为
\begin{eqnarray}
a'_{pp} &=& a_{pp} - t a_{pq}, \label{jacobi-a'pp-2}, \\
a'_{qq} &=& a_{qq} + t a_{pq}, \label{jacobi-a'qq-2}, \\
a'_{ip} &=& a_{ip} - s (a_{iq} + \frac{s}{1+c} a_{ip}), \\
a'_{iq} &=& a_{iq} + s (a_{ip} - \frac{s}{1+c} a_{iq}).
\end{eqnarray}

到这里，我们理清了旋转变换怎样使$a'_{pq}$等于零。
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
我们可以进行一系列这样的旋转变换。
有时变为零的非对角元，可能在下一次旋转变换之后，又变为非零。
比如
\begin{equation}
A'' = P^\top_{pr} P^\top_{pq} A P_{pq} P_{pr},
\end{equation}
经过第一次正交变换$P_{pq}$，$a'_{pq}=0$，但是经过第二次正交变换$P_{pr}$，根据\\
{\centering \color{blue}(\ref{jacobi-a'ip}): $a'_{ip} = c a_{ip} - s a_{iq}, ~~\forall i \neq p,q $ }\\
有
\begin{equation}
a''_{pq} = a''_{qp} = ca'_{qp} - sa'_{qq} = -s a'_{qq} \neq 0.
\end{equation}
那么这一系列正交变换会使得A最终变为对角阵吗？
答案是肯定的。
\end{frame}

\begin{frame}[fragile]{Jacobi方法：正交变换后，矩阵所有矩阵元平方和不变。}
这是正交变换的性质，这很容易证明。
矩阵$N=A^\top A$的对角元为
\begin{equation}
N_{ii} = \sum_j A^\top_{ij} A_{ji} = \sum_j A_{ji} A_{ji},
\end{equation}
所以$N$的迹为
\begin{equation}
Trace(N) = \sum_i N_{ii} = \sum_{ij} A_{ji} A_{ji},
\end{equation}
正是矩阵$A$的所有矩阵元的平方和。
经过正交变换$P$以后，矩阵$A$变为
\begin{equation}
A'=P^\top A P,
\end{equation}
它的所有矩阵元平方和为
\begin{eqnarray}
Trace( A'^\top A) &=& Trace( P^\top A^\top P P^\top A P ) = Trace( P^\top A^\top A P ) 
\nonumber\\
&=& Trace( P^\top N P) = Trace(N) = \sum_{ij} A_{ji} A_{ji}.
\end{eqnarray}
所以正交变换前后，矩阵的所有矩阵元平方和不变。
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
根据(\ref{jacobi-a'pp}-\ref{jacobi-a'iq})，正交变换$P_{pq}$之后，
\begin{equation}
a'^2_{pp} + a'^2_{qq} = a^2_{pp} + a^2_{pq} + a^2_{qp} + a^2_{qq},
\end{equation}
所以矩阵$A$经过这次变换以后，对角元平方和增加了 $2 a^2_{pq}$。

对角元平方和的上限为$A$的所有矩阵元平方和，而我们可以任意选取$p,q$进行$P_{pq}$变换，所以最终一定可以将$A$对角化。
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
如果把所有相似变换矩阵$P_{pq}$的乘积记作
\begin{equation}
S = P_{p_1 q_1} P_{p_2 q_2} \cdots P_{p_N q_N},
\end{equation}
对角化后的矩阵记为$\Lambda$，则有
\begin{equation}
S^\top A S = \Lambda,
\end{equation}
即
\begin{equation}
AS = S \Lambda,
\end{equation}
也就是说，S的各列是A的特征向量，$\Lambda$的对角元是A的特征值。
所以在程序中，可以一边做正交变换，一边记下它们的总乘积，即同时得到特征值和特征向量。
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
实现的代码如下
\begin{lstlisting}
/*
 * jacobi(...): 用Jacobi算法对角化一个对称实矩阵，得到本征值和本征矢
 *
 * int n        矩阵维数
 * double *a    实对称阵A，a[i*n+j]表示A_{ij}，i,j=0,1,...,n-1
 * double *z    本征值：z[0],z[1],...,z[n-1]
 * 程序运行结束以后，本征矢量储存在 a 中，覆盖原来的输入矩阵，a的第i行是对应于z[i]的本征矢量，即a[ i*n +j ], j=0,...,n-1 是第i个本征矢量。
 */
void jacobi(int n, double *a, double *z)
{
    int l, i, j, p, q;
    double max, temp, s, c, theta, t;
    double g, h;

    double **v = new double *[n];//v矩阵初始时是一个nxn的单位阵，经过for(l)循环以后，v是所有相似矩阵的乘积，其列向量构成原来的矩阵 a 的特征向量。
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
\lstset{firstnumber = last}
\begin{lstlisting}
    for(i=0;i<n;i++){
        v[i] = new double [n];
        for(j=0;j<n;j++){
            if(j!=i) v[i][j]=0;
            else v[i][j]=1;
        }
    }

    for (l = 0; l < max_iteration; l++){
        max = 0;// 找到绝对值最大的非对角元,记下它的行数p、列数q
        for (i = 0; i < n; i++){
            for (j = 0; j < n; j++){
                if ( i != j && fabs( a[i*n+j] ) >max ){
                    max = fabs( a[i*n+j] );
                    p = i;
                    q = j;
                }
            }
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
\lstset{firstnumber = last}
\begin{lstlisting}
        if ( max < 1e-10 ) break;
        //如果最大的非对角元已经小于1E-10，说明矩阵已经对角化了，所以跳出for(l)循环

        //构造相似变换
        theta = ( a[q*n+q] - a[p*n+p] )/(2*a[p*n+q]); //参考课件(70)式计算theta
        t = 1/ ( fabs(theta) + sqrt(theta*theta+1) );// 参考课件(75)式计算t
        if(theta<0) t *= -1;

        c = 1/sqrt(1+t*t); //参考课件(76)式计算c,s
        s = c*t;

        a[p*n+p] -= t* a[p*n+q];//参考课件(77-80)计算a'pp, a'qq, a'ip, a'iq
        a[q*n+q] += t* a[p*n+q];
        a[p*n+q] = 0;
        a[q*n+p] = 0;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
\lstset{firstnumber = last}
\begin{lstlisting}
        for(i=0;i<n;i++){//计算第p,q行，或第p,q列其它元素
            if(i!=p&&i!=q){
                temp = a[i*n+p];
                a[i*n+p] -= s*(a[i*n+q] + s/(1+c)*temp);
                a[i*n+q] += s*(temp - s/(1+c)*a[i*n+q]);
                a[p*n+i] = a[i*n+p];
                a[q*n+i] = a[i*n+q];
            }
        }
        for (i = 0; i < n; i++){  // v = P1 P2 ... 
            temp = v[i][p];
            v[i][p] =  temp * c - v[i][q] * s;
            v[i][q] =  temp * s + v[i][q] * c; 
        }
        /* for test only
        cout<<"a:"<<endl;
        for(i=0;i<n;i++){
            for(j=0;j<n;j++)
                cout<<a[i*n+j]<<",";
            cout<<endl;
        }
        */
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
\lstset{firstnumber = last}
\begin{lstlisting}
    }//l
    if( l== max_iteration){//如果经过for(l)循环以后，l的值为max_iteration，说明矩阵仍然没有对角化，我们就报告没有收敛，然后退出程序
        cout<<" Jacobi: after "<<max_iteration<<" iterations, didn't converge...\n";
        exit(1);
    }
    for (i = 0; i < n; i++)z[i] = a[i*n+i];//记录本征值到z中
    for(i=0;i<n;i++){ //把矩阵v中的数据录到a中，录完以后，a[i*n+0]~a[i*n+n-1]等于v的第i列
        for(j=0;j<n;j++){
            a[i*n+j] = v[j][i];
        }
    }
    for(i=0;i<n;i++) delete [] v[i]; delete [] v;//释放v的动态内存
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
我们可以写一个主函数，调用jacobi函数，测试这个函数，验证它是否正确。
\lstset{firstnumber = last}
\begin{lstlisting}
#include<iostream>
using namespace std;

#include<cmath>
#include<ctime>

#include"library.h"

int main(){

    srand(time(NULL));
    int i,j,k;
    int n=3;
    double y;

    double * a = new double [n*n];
    double *a_diag = new double [n*n];
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
\lstset{firstnumber = last}
\begin{lstlisting}
    for(i=0;i<n;i++){// 给a和a_diag赋值，a_diag作为参数传入jacobi(...)
        for(j=i;j<n;j++){
            a[ i*n + j ] = (i+2)*(j+2)+3;//(double)rand()/RAND_MAX;
            a[ j*n + i ] = a[ i*n + j];
            a_diag[ i*n + j] = a[ i*n + j];// a_diag = a
            a_diag[ j*n + i] = a[ j*n + i];
        }
    }
    cout<<"A:"<<endl;// 输出最初的矩阵A
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            a_diag[ i*n+j ] = a[ i*n+j ];
            cout<<a_diag[ i*n+j ]<<",";
        }
        cout<<endl;
    }

    double *z = new double [n];

    jacobi(n, a_diag, z);//调用jacobi函数，将a_diag对角化
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
\lstset{firstnumber = last}
\begin{lstlisting}
    //输出本征值
    cout<<"jacobi eigenvalues: "; for(int i=0;i<n;i++)cout<<z[i]<<","; cout<<endl;
    //输出 A \vec{v}，即用原来的矩阵A乘以求得的本征矢量
    cout<<"A v: "<<endl;
    for(i=0;i<n;i++){// ith column of ( A (a_diag)^top )
        cout<<" A v"<<i<<" = ";
        for(j=0;j<n;j++){
            y=0;
            for(k=0;k<n;k++) // A \vec{v} 
                y += a[j*n+k] * a_diag[i*n+k];
            cout<<y<<",";
        }
        cout<<endl;
    }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Jacobi方法}
\lstset{firstnumber = last}
\begin{lstlisting}
    //输出 lambda * \vec{v}，即用求得的本征值乘以本征矢量
    //如果计算正确，就会得到 A \vec{v} = \lambda \vec{v}，否则程序就写错了。
    cout<<" lambda * v : "<<endl;
    for(i=0;i<n;i++){
        cout<<" lambda"<<i<<" * v"<<i<<" = ";
        for(j=0;j<n;j++)cout<< z[i]* a_diag[ i*n + j]<<",";
        cout<<endl;
    }

    delete [] z,a,a_diag;
    return 0;
}
\end{lstlisting}
\end{frame}
\begin{frame}[fragile]{Householder变换}
Householder变换是一种正交变换，通过它可以将对称阵转化为三对角矩阵，三对角矩阵就容易处理了。

Householder变换由一系列正交相似变换构成,
\begin{eqnarray}
S &=& S_1 S_2 \cdots S_{n-2}, \\
S_1 S^\top_1 &=& S_2 S^\top_2 = \cdots = 1, \\
S^\top A S &=& 
\begin{pmatrix}
a_{11} & e_2 &    &  & \\
e_2 & a'_{22} & e_2 &  & \\
    & \ddots & \ddots & \ddots \\
    &     &         & a^{(n-2)}_{n-1,n-1} & e_{n-1} \\
    &     &         & e_{n-1}             & a^{(n-2)}_{n,n}
\end{pmatrix}
\end{eqnarray}
上面的$a'_{22}$经过了一次迭代（一次相似变换），$a^{(n-2)}_{n-1,n-1}, a^{(n-2)}_{n,n}$经过了$n-2$次迭代。
\end{frame}

\begin{frame}[fragile]{Householder变换}
$S_1$这样构造：
\begin{eqnarray}
S_1 = \begin{pmatrix}
1 & \vec{0}^\top \\
\vec{0} & P
\end{pmatrix},
\end{eqnarray}
其中$\vec{0}$表示n-1个0构成的列矢量，$\vec{0}^\top$是它的转置，$P$满足$P^\top = P, P^2 = 1$。
这样构造：
\begin{equation}
P = 1 - 2 \vec{u} \vec{u}^\top,
\end{equation}
其中$\vec{u}$是归一化的矢量，那么P的这个形式就可以保证$P^\top = P, P^2 = 1$。

为了方便叙述，我们把A分为四块矩阵
\begin{equation}
A = \begin{pmatrix}
a_{11} & \vec{v}^\top \\
\vec{v} & B
\end{pmatrix},
\end{equation}
其中，$\vec{v} = \{ a_{21}, \cdots, a_{n1} \}^\top$。
\end{frame}

\begin{frame}[fragile]{Householder变换}
经过第一次相似变换以后，$A$变为
\begin{equation}
S^\top_1 A S_1 =
\begin{pmatrix}
a_{11} & (P \vec{v})^\top \\
P\vec{v} & B'
\end{pmatrix},
\end{equation}
其中$B' = PAP$。

所以我们想，如果$P\vec{v}$只有第一个元素不为零，其它都是零，我们就离目标（三对角阵）接近了一步。

$(P \vec{v})$的模平方为
\begin{equation}
(P \vec{v})^\top (P\vec{v}) = |\vec{v}|^2,
\end{equation}
所以，我们得命令
\begin{equation}
P \vec{v} = \pm v \vec{e}_2,
\end{equation}
其中v是$\vec{v}$的模，$\vec{e}_2$是$\vec{v}$的线性空间第一个单位矢量。
\end{frame}

\begin{frame}[fragile]{Householder变换}
这样，我们就把$\vec{u}$定下来了，
\begin{equation}
P\vec{v} = (1- 2 \vec{u} \vec{u}^\top) \vec{v}
= \vec{v} - 2 \vec{u} (\vec{u}^\top \vec{v}) 
= \pm v \vec{e}_2,
\end{equation}
最后一个等号两边同时与$\vec{v}$做点乘，得到
\begin{equation}
2(\vec{u}^\top \vec{v} )^2 = v^2 \mp a_{21} v.
\end{equation}
所以可以定义
\begin{equation}
\vec{u} = \frac{ \vec{v} \mp v \vec{e}_2 }{ 2 \vec{u}^\top \vec{v} } 
= \frac{\vec{v} \mp v \vec{e}_2}{\sqrt{ 2v^2 \mp 2a_{21} v  }},
\end{equation}
其中$\mp$对应两种选择，都可以完成目标，但我们一般选$v^2 \mp a_{21} v$更大的那一种，防止误差太大。
\end{frame}

\begin{frame}[fragile]{Householder变换}
\begin{shaded*}
这里我们可以做一个观察，如果把$\vec{u}$定为
\begin{equation}
\vec{u} = \frac{\vec{v} \mp v \vec{e}_k}{\sqrt{ 2v^2 \mp 2a_{k1} v  }}, ~~~ k=2,3,...,n
\end{equation}
那么，$P=1-2\vec{u}\vec{u}^\top$作用到$\vec{v} = \{ a_{21}, a_{31}, \cdots, a_{n1} \}^\top$上，会得到
\begin{equation}
P\vec{v} = \vec{v} - 2 \vec{u} (\vec{u}^\top \vec{v}) = \pm v \vec{e}_k.
\end{equation}
也就是，$S^\top_1 A S_1$第一列只有第1个和第k个元素不为零，第一行只有第1个和第k个元素不为零。
所以，这个技巧可以针对$k=2,...,n$的，使得$\vec{v}$中的相应的分量不为零，其余变为零。
这样对一个矢量$\vec{v}$做变换，原本就叫做Householder变换，它实际上是将矢量$\vec{v}$相对一个超平面做镜像反射变换，这个超平面的法向矢量为$\vec{u}$。
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{Householder变换}
\begin{shaded*}
我们可以再做一个观察，如果把$S_1$定为
\begin{eqnarray}
S_1 = \begin{pmatrix}
P & \vec{0}^\top \\
\vec{0} & 1
\end{pmatrix},
\end{eqnarray}
矩阵A划分为
\begin{equation}
A = \begin{pmatrix}
B & \vec{v} \\
\vec{v}^\top & a_{nn}
\end{pmatrix},
\end{equation}
会发现
\begin{equation}
\begin{pmatrix}
A' & P \vec{v}\\
(P\vec{v})^\top & a_{nn}
\end{pmatrix},
\end{equation}
其中$\vec{v} = \{ a_{1n}, a_{2n}, \cdots, a_{n-1,n} \}$。
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{Householder变换}
\begin{shaded*}
然后，我们可以定义
\begin{equation}
\vec{u} = \frac{\vec{v} \mp v \vec{e}_{n-1}}{\sqrt{ 2v^2 \mp 2a_{n-1,1} v  }}, ~~~ k=2,3,...,n
\end{equation}
这样，$S^\top_1 A S_1$的最后一行、最后一列就只有第$n-1,n$个元素不为零，其它都是零。

所以，我们这一节介绍的Householder变换，也可以从最后一行（列）开始，逐渐向三对角形式靠近。
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{Householder变换}
把定下的$\vec{u}$带入P，就得到了$S_1$，然后进行相似变换，就完成了第一步。
但是直接做两次矩阵乘法，太费时间了，我们可以如下简化计算
\begin{eqnarray}
\vec{p} &=& 2 B \vec{u}, \\
BP &=& B(1-2\vec{u}\vec{u}^\top) = B - \vec{p} \vec{u}^\top, \\
PBP &=& (1-2\vec{u}\vec{u}^\top)(B - \vec{p} \vec{u}^\top) \nonumber\\
    &=& B - \vec{u}\vec{p}^\top - \vec{p}\vec{u}^\top + 2 (\vec{u}^\top \vec{p}) \vec{u}\vec{u}^\top.
\end{eqnarray}
如果定义
\begin{eqnarray}
K &\equiv& \vec{u}^\top \vec{p}, \nonumber\\
\vec{q} &=& \vec{p} - K \vec{u},
\end{eqnarray}
则有
\begin{eqnarray}
PBP = B - \vec{u}\vec{q}^\top - \vec{q}\vec{u}^\top.
\end{eqnarray}

说到这里，我们解释了$S_1$的构造思路，并且展示了效果以及计算方法。
$S_2$则只作用在$B'$上，将$B'$第一行、第一列第$3,\cdots,n-1$个元素变为零。
如此继续下去，通过$n-2$次相似变换，最终将$A$转换为三对角阵。
\end{frame}

\begin{frame}[fragile]{一维无限深势井}
在第八章线性方程组的开头，我们提到过，利用有限差分，二阶常微分方程可以转化为线性方程组求解问题。
那么，如果微分方程中有待定的特征值，则可以考虑将微分方程转化为矩阵本征值问题。
1维无限深势井，就是一个这样的例子。
\begin{equation}
V(x) = \left\{
\begin{aligned}
&0, ~~~ 0<x<a\\
&\infty, ~~else
\end{aligned}
\right.
\end{equation}
在$[0,a]$区间内，薛定谔方程为
\begin{equation}
- \frac{ \hbar^2 }{2m} \frac{ d^2 }{dx^2} \psi = E \psi，
\end{equation}
相应的本征值和解为
\begin{eqnarray}
E_n &=& \frac{ n^2 \pi^2 \hbar^2 }{ 2 m a^2}, \\
\psi_n &=& \sqrt{\frac{2}{a}} \sin (n \pi x /a  ).
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{一维无限深势井}
在数值中，为了简便，我们取$\frac{\hbar^2}{2m} = a = 1$，则薛定谔方程变为
\begin{equation}
- \frac{d^2}{dx^2} \psi = E \psi, \label{linear-1d-infinite-well}
\end{equation}
其解为
\begin{eqnarray}
E_n &=& n^2 \pi^2 , \\
\psi_n &=& \sqrt{2} \sin (n \pi x  ).
\end{eqnarray}
我们可以将$[0,a]$，即$[0,1]$区间等分$n+1$份，则有$n$个节点$x_1, \cdots, x_{n}$，当然，还有$x_{0} = 0, x_{n+1} = 1$。
相应地，记各个节点上$\psi$的数值近似为$y_{0}, y_1, ..., y_{n+1}$，由边界条件，我们有
\begin{equation}
y_{0} = 0, ~~ y_{n+1} = 0.
\end{equation}
\end{frame}

\begin{frame}[fragile]{一维无限深势井}
利用有限差分公式
\begin{eqnarray}
\frac{d^2}{dx^2} \psi (x) |_{x_i} &\approx& ( \psi (x_{i+1}) + \psi(x_{i-1}) - 2\psi(x_i) )/h^2, 
\end{eqnarray}
得到$y_0, \cdots, y_{n+1}$之间的关系。
在$x_1, x_2, \cdots, x_n$各个节点上，(\ref{linear-1d-infinite-well})转化为
\begin{eqnarray}
-\frac{1}{h^2}y_0 + \frac{2}{h^2}y_1 - \frac{1}{h^2}y_2 &=& E y_1, \\
-\frac{1}{h^2}y_1 + \frac{2}{h^2}y_2 - \frac{1}{h^2}y_3 &=& E y_2, \\
-\frac{1}{h^2}y_2 + \frac{2}{h^2}y_3 - \frac{1}{h^2}y_4 &=& E y_3, \\
\cdots&& \\
-\frac{1}{h^2}y_{n-1} + \frac{2}{h^2}y_n - \frac{1}{h^2}y_{n+1} &=& E y_n, 
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{一维无限深势井}
将$y_0 = y_{n+1} = 0$代入其中，就得到一个矩阵形式的本征值问题：
\begin{equation}
\begin{bmatrix}
\frac{2}{h^2}  & - \frac{1}{h^2}  \\
-\frac{1}{h^2} & \frac{2}{h^2}   & - \frac{1}{h^2}\\
				& -\frac{1}{h^2} & \frac{2}{h^2}   & - \frac{1}{h^2} \\
				&				 & \vdots          & \vdots        & \vdots \\
				&				&				&	- \frac{1}{h^2} &  	\frac{2}{h^2}	
\end{bmatrix}
\begin{bmatrix}
y_1 \\
y_2 \\
y_3 \\
\vdots \\
y_n
\end{bmatrix}
= E
\begin{bmatrix}
y_1 \\
y_2 \\
y_3 \\
\vdots \\
y_n
\end{bmatrix}
\end{equation}
对角化上式左侧的矩阵，就得到了能量和波函数。
\end{frame}

\begin{frame}[fragile]{一维无限深势井}
下面的函数，即实现上面的算法，
\begin{lstlisting}
/* ISW1d:       解出1维无限深方势井本征值和波函数。假定 hbar^2/2m = a = 1，ISW1d: Infinite-Square-Well-1-dimension。解出的波函数写进"ISW1d_wave.txt"中。
 * int n:       (0,1)上的节点数
 * double *e:   本征值，从小到大排列
 */ 
double ISW1d(int n, double *e){//1d Infinite Square Well
        int i,j;
        double h = 1.0/(n+1);// 将区域 [0,1] 分为 n+1 份，所以其中有 n 个节点, 步长为h
        double *a = new double [n*n];//课件中的矩阵
        for(i=0;i<n;i++){
                for(j=0;j<n;j++){
                        if(i==j){
                                a[i*n+j] = 2/h/h;
                        }
                        else if( abs(i-j)==1 ) {
                                a[i*n+j] = -1/h/h;
                        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{一维无限深势井} 
\lstset{firstnumber = last}
\begin{lstlisting}
                        else a[i*n+j]=0;
                }
        }       
        jacobi(n, a, e);//对角化
        sort_eig(n, a, e);//按本征值从小到大，将本征值和本征矢量排序
        // 将波函数输出到ISW1d_wave.txt
        ofstream fp("ISW1d_wave.txt");
        for(i=0;i<n+1;i++)fp<<"0.0\t"; fp<<endl;
        for(i=0;i<n;i++){
                fp<<setprecision(4)<<(i+1)*h<<"\t";
                for(j=0;j<n;j++)
                        fp<<setprecision(4)<<a[j*n+i]<<"\t";
                fp<<endl;
        }
        fp<<"1.0\t";
        for(i=0;i<n;i++)fp<<"0.0\t"; fp<<endl;
        fp.close();

        delete [] a;//释放动态内存        
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{一维无限深势井} 
这个函数可以求出本征值来。
如果要得到更精确的结果，我们需要增大$n$，减小步长$h$。
但是计算量随着$n$的变大，迅速地变大，是否有更好的办法？

我们可以这样看待这个问题。
求得的$E_i$都是步长$h$的函数，当$h$趋于0的时候，$E_i$一定会趋于理论值，只是我们算不动太大的矩阵。
那么，我们可以选取不同的步长值，$h, h/2, h/4$，分别得到$E^h_i, E^{h/2}_i, E^{h/4}_i$，然后做多项式外推，推测$h=0$时$E_i$的值。
\end{frame}

\begin{frame}[fragile]{一维无限深势井} 
在下面的代码中，我们调用第二部分写的Nevill函数，用Nevill外推法推得$E_i$的值。
\begin{lstlisting}
int main(){
        
        int n=30, i, j;
        double y, dy;

        double x_data[3], eig[10][3];//用作Nevill外推的输入量

        for(i=0;i<3;i++){
                x_data[i] = 1.0/(n+1);//步长
                double *e = new double [n];
                ISW1d(n, e);
                cout<<"e:"; for(j=0;j<10;j++)cout<<e[j]/M_PI/M_PI<<","; cout<<endl;
                for(j=0;j<10;j++) eig[j][i] = e[j];//本征值
                delete [] e;
                n*=2;
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{一维无限深势井} 
\lstset{firstnumber = last}
\begin{lstlisting} 
        cout<<"Extrapolation:"<<endl;
        for(i=0;i<10;i++){
                Nevill(3, x_data, eig[i], 0.0, y, dy);// y为外推得到的 E_{i+1}
                cout<<"E"<<(i+1)<<"/pi/pi = "<<y/M_PI/M_PI<<endl;
        }
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{一维无限深势井} 
外推得到的前10个能级如下
\begin{table}[h]
\centering
\begin{tabular}{ccc}
\hline\hline
i & $E_i/\pi^2$ & 理论值 \\
\hline
1 & 1.000		& 1 \\
2 & 4.000		& 4 \\
3 & 9.000		& 9 \\
4 & 16.000		& 16 \\
5 & 25.001		& 25 \\
6 & 36.003		& 36 \\
7 & 49.008		& 49 \\
8 & 64.017		& 64 \\
9 & 81.035		& 81 \\
10 & 100.066 	& 100 \\
\hline\hline
\end{tabular}
\caption{一维无限深势阱前10个能级，保留三位小数。利用有限差分，将薛定谔方恒转化为矩阵本征值问题，然后用不同的步长计算$E_i$，再外推得到步长$h\rightarrow 0$时的$E_i$值。}
\end{table}
\end{frame}

\begin{frame}[fragile]{一维无限深势井} 
上面的计算中，用到的矩阵的最大维度为$120 \times 120$，但已经推得非常精确的结果。
下图中是$n=120$时求得的最低的几个能级的波函数
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{ISW1d}
%\caption{一维无限深势井的波函数。}
\end{figure}
\end{frame}

\begin{frame}[fragile]{作业}
1. 怎样用高斯消元法求矩阵的逆矩阵？是否可以求本征值？（简答，将回答写在邮件正文即可）\\
2. 利用ludcmp函数，求矩阵的行列式：
\begin{equation}
A=
\begin{pmatrix}
1   & 1 &   1   &   1\\
1   & 2 &   2   &   2\\
1   & 2 &   3   &   3\\
1   & 2 &   3   &   4
\end{pmatrix}.
\end{equation}
*3. 利用ludcmp函数，求矩阵的逆矩阵：
\begin{equation}
A=
\begin{pmatrix}
1   & 1 &   1   &   1\\
1   & 2 &   2   &   2\\
1   & 2 &   3   &   3\\
1   & 2 &   3   &   4
\end{pmatrix}.
\end{equation}
4. 实验室制备了100个完全相同的放射性X原子，它们在单位时间内衰变的概率相同。每隔十分钟数一遍剩下的 X 原子数,得到数据：
(9:00, 100), (9:10, 71), (9:20, 50), (9:30, 37), 
(9:40, 26), (9:50, 17), (10:00, 11), (10:10, 9)
试用最小二乘法拟合这种原子的半衰期。
\end{frame}

\begin{frame}[fragile]{作业}
**5. 试编写jacobi迭代法程序，求解以下方程组的解
\begin{equation}
\begin{pmatrix}
5   & 1 &   1   &   1\\
1   & 6 &   2   &   2\\
1   & 2 &   8   &   3\\
1   & 2 &   3   &   10
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3 \\
x_4
\end{pmatrix}
=\begin{pmatrix}
8\\
11\\
14\\
16
\end{pmatrix}.
\end{equation}

6. 利用jacobi方法，求以下矩阵的本征值与本征矢量：
\begin{equation}
A = \begin{pmatrix}
0.575155 & 0.878075 & 0.939033\\
0.878075 & 0.445565 & 0.99726\\
0.939033 & 0.99726 & 0.957276
\end{pmatrix}
\end{equation}
% 答案：eigenvalues: -0.397167,-0.183899,2.55906,
%eigenkets: 
%-0.482712,0.830803,-0.277047,
%0.687128,0.163144,-0.707983,
%0.542996,0.532118,0.649619,
\end{frame}

\begin{frame}[fragile]{作业}
7. 试用课件中展示的有限差分+矩阵对角化+Nevill外推方法，求解1维谐振子势
\begin{equation}
V(x) = \frac{1}{2} k x^2,
\end{equation}
为了叙述简便，不妨取$\hbar = m = k = 1$, 则薛定谔方程为
\begin{eqnarray}
- \frac{d^2}{dx^2} \psi + x^2 \psi &=& 2 E \psi, \\
                                 E &=& (n+\frac{1}{2}).
\end{eqnarray}
我们可以选取区间$[-5,5]$，将区间等分为$n+1$份，并假设$\psi(-5) = \psi(5) = 0$，然后进行类似于讲义中类似于一维无限深势井的处理，\\
1)试推得前3个本征值，\\
*2)画出前三个本征态的波函数，并做图。
\end{frame}

\end{document}